#ifndef AMCUTILITYALGFFT_H
#define AMCUTILITYALGFFT_H

#include <math.h>

/* 自定义复数double类型结构体 */
typedef struct TComplexNumber_struct{
    double dReal;    //实部
    double dImag;	 //虚部
}TComplexNumber;

class AMCUtilityAlgFFT
{
public:

	AMCUtilityAlgFFT() {};

    /******************************************************************************
    * 模块名称：real_fft_power_dB(double pdInputData[], int iInputDataLength, double pdOutPutFftPower[], int *iOutputBufferLength, int iWantedFftSize)
    * 模块功能：计算输入实数序列的功率谱取dB
    * 参数说明： pdInputData[]          - double类型输入数据buffer指针
    *		    iInputDataLength       - 输入数据长度
    *		    pdOutPutFftPower[]     - double类型输出功率谱buffer指针
    *		    *iOutputBufferLength   - 输出功率谱buffer长度，计算后修改该参数为buffer中有效数据点数
    *		    iWantedFftSize         - 输入希望的FFT计算点数，若输入数据不够，末端补0
    * 返回值  ：执行成功标识，返回0表示执行成功，返回-1表示实际生成的功率谱点数超过了功率谱buffer深度，错误
    * 建立时间：2020-08-12
    * 创建人员：陈文棵
    * 修改日期：
    * 修改人员：
    *******************************************************************************/
    template< class T>
    static int real_fft_power_dB(T pdInputData[], int iInputDataLength, T pdOutPutFftPower[], int *iOutputBufferLength, int iWantedFftSize);
};

#define Pi 3.14159265358979323846264338

/******************************************************************************
* 模块名称：tcn_mul(TComplexNumber tCNa, TComplexNumber tCNb, TComplexNumber *tCNc)
* 模块功能：复数乘法
* 参数说明：tCNa   - TComplexNumber复数类型乘数1
*		    tCNb   - TComplexNumber复数类型乘数2
*	        *ptCNc - TComplexNumber复数指针类型乘积
* 返回值  ：执行成功标识，返回0表示执行成功
* 建立时间：2020-08-12
* 创建人员：陈文棵
* 修改日期：
* 修改人员：
*******************************************************************************/
static void tcn_mul(TComplexNumber tCNa, TComplexNumber tCNb, TComplexNumber *ptCNc)
{
    ptCNc->dReal = tCNa.dReal * tCNb.dReal - tCNa.dImag * tCNb.dImag;
    ptCNc->dImag = tCNa.dReal * tCNb.dImag + tCNa.dImag * tCNb.dReal;
}

/******************************************************************************
* 模块名称：tcn_sub(TComplexNumber tCNa, TComplexNumber tCNb, TComplexNumber *tCNc)
* 模块功能：复数减法
* 参数说明：tCNa   - TComplexNumber复数类型被减数
*		    tCNb   - TComplexNumber复数类型减数
*	        *ptCNc - TComplexNumber复数指针类型差
* 返回值  ：执行成功标识，返回0表示执行成功
* 建立时间：2020-08-12
* 创建人员：陈文棵
* 修改日期：
* 修改人员：
*******************************************************************************/
static void tcn_sub(TComplexNumber tCNa, TComplexNumber tCNb, TComplexNumber *ptCNc)
{
    ptCNc->dReal = tCNa.dReal - tCNb.dReal;
    ptCNc->dImag = tCNa.dImag - tCNb.dImag;
}

/******************************************************************************
* 模块名称：tcn_plus(TComplexNumber tCNa, TComplexNumber tCNb, TComplexNumber *tCNc)
* 模块功能：复数加法
* 参数说明：tCNa   - TComplexNumber复数类型加数1
*		    tCNb   - TComplexNumber复数类型加数2
*	        *ptCNc - TComplexNumber复数指针类型和
* 返回值  ：执行成功标识，返回0表示执行成功
* 建立时间：2020-08-12
* 创建人员：陈文棵
* 修改日期：
* 修改人员：
*******************************************************************************/
static void tcn_plus(TComplexNumber tCNa, TComplexNumber tCNb, TComplexNumber *ptCNc)
{
    ptCNc->dReal = tCNa.dReal + tCNb.dReal;
    ptCNc->dImag = tCNa.dImag + tCNb.dImag;
}

/******************************************************************************
* 模块名称：get_log2N(int iN, int *iLog2N, int *iDvalue)
* 模块功能：对于正整数iN，计算log2N，当iN不是2的整数次方时，iLog2N为大于log2N且最接近的整数，iDvalue为2^iLog2N - iN的值
* 参数说明：iN       - 计算对象
*		    *iLog2N  - 大于等于log2N且最接近的整数
*	        *iDvalue - 2^iLog2N - iN
* 返回值  ：执行成功标识，返回0表示执行成功, 返回-1表示iN不是正整数
* 建立时间：2020-08-12
* 创建人员：陈文棵
* 修改日期：
* 修改人员：
*******************************************************************************/
static int get_log2N(int iN, int *iLog2N, int *iDvalue)
{
    int i;
    if (iN <= 0)
    {
        return -1;
    }

    for (i = 30; i >= 0; i--)
    {
        if (iN > (1 << i))
        {
            *iLog2N = i + 1;
            *iDvalue = (1 << (*iLog2N)) - iN;
            break;
        } /* end of for if (N > (1 << i)) */
        else if (iN == (1 << i))
        {
            *iLog2N = i;
            *iDvalue = 0;
            break;
        } /* end of for else if (N == (1 << i)) */
    } /* end of for for (int i = 30; i >= 0; i--) */

    return 0;
}

/******************************************************************************
* 模块名称：get_fftinput_reshape_idx(int iFftSize, int iFftSteps, int piFftInputReshapeIdx[])
* 模块功能：计算对FFT的输入数据进行重新排序的索引隐射表
* 参数说明：iFftSize  - FFT点数
*		    iFftSteps - FFT蝶形运算的级数
*	        piFftInputReshapeIdx[] - 重新排序的索引buffer指针
* 返回值  ：执行成功标识，返回0表示执行成功
* 建立时间：2020-08-12
* 创建人员：陈文棵
* 修改日期：
* 修改人员：
*******************************************************************************/
static int get_fftinput_reshape_idx(int iFftSize, int iFftSteps, int piFftInputReshapeIdx[])
{
    int iTmpIdx = 0;
    int i,j;
    for (i = 0; i < iFftSize; i++)
    {
        iTmpIdx = 0;
        for (j = 0; j < iFftSteps; j++)
        {
            iTmpIdx = iTmpIdx + (((i >> j) & 1) << (iFftSteps - j - 1));
        } /* end of for for for (int j = 0; j < N_fft_steps; j++) */
        piFftInputReshapeIdx[i] = iTmpIdx;
    } /* end of for for (int i = 0; i < fft_size; i++) */

    return 0;
}

/******************************************************************************
* 模块名称：get_fft_wn(int iFftSize, TComplexNumber ptCNFftWn[])
* 模块功能：计算FFT各级蝶形运算需要使用的旋转因子
* 参数说明：iFftSize    - FFT点数
*		    ptCNFftWn[] - 存储旋转因子buffer的指针，buffer深度为FFT点数一半
* 返回值  ：执行成功标识，返回0表示执行成功
* 建立时间：2020-08-12
* 创建人员：陈文棵
* 修改日期：
* 修改人员：
*******************************************************************************/
static int get_fft_wn(int iFftSize, TComplexNumber ptCNFftWn[])
{
    int i;
    for (i = 0; i < iFftSize / 2; i++)
    {
        ptCNFftWn[i].dReal = cos(2 * Pi*i / iFftSize);
        ptCNFftWn[i].dImag = -sin(2 * Pi*i / iFftSize);
    } /* end of for for (int i = 0; i < iFftSize / 2; i++) */

    return 0;
}

/******************************************************************************
* 模块名称：fft_butterfly(TComplexNumber ptCNInputData[], TComplexNumber ptCNOutputData[], TComplexNumber ptCNFftWn[], int piReshapeIdx[], int iFftSize, int iFftSteps)
* 模块功能：FFT算法核心--多级蝶形运算
* 参数说明：ptCNInputData[]  - TComplexNumber复数类型蝶形运算第一级输入数据序列buffer指针
*		    ptCNOutputData[] - TComplexNumber复数类型蝶形运算最后一级输出数据序列buffer指针
*		    ptCNFftWn[]      - TComplexNumber复数类型蝶形运算需要使用的输入旋转因子序列buffer指针
*		    piReshapeIdx[]   - int类型蝶形运算需要使用的输入数据排序索引序列buffer指针
*		    iFftSize         - int类型输入FFT总点数
*		    iFftSize         - int类型输入FFT蝶形运算级数
* 返回值  ：执行成功标识，返回0表示执行成功
* 建立时间：2020-08-12
* 创建人员：陈文棵
* 修改日期：
* 修改人员：
*******************************************************************************/
static int fft_butterfly(TComplexNumber ptCNInputData[], TComplexNumber ptCNOutputData[], TComplexNumber ptCNFftWn[], int piReshapeIdx[], int iFftSize, int iFftSteps)
{
    int iFftPointNumPerBufferfly = 0;
    int iFftButterflyNumPerStep  = 0;
    int iTempStartIdx = 0;
    int i,j,k;
    TComplexNumber tCNTempWn;
    TComplexNumber tCNTempMultData;
    TComplexNumber *ptCNTempFftInputData;
    TComplexNumber *ptCNTempFftOutputData;
    TComplexNumber *ptCNTempPointer;
    ptCNTempFftInputData = new TComplexNumber[iFftSize];
    //ptCNTempFftInputData = malloc(iFftSize * sizeof(TComplexNumber));
    ptCNTempFftOutputData = new TComplexNumber[iFftSize];
    //ptCNTempFftOutputData = malloc(iFftSize * sizeof(TComplexNumber));

    for (i = 0; i < iFftSize; i++)
    {
        ptCNTempFftInputData[i] = ptCNInputData[piReshapeIdx[i]];
    } /* end of for for (int i = 0; i < iFftSize; i++) */
    for (i = 0; i < iFftSteps; i++)
    {
        iFftPointNumPerBufferfly = (1 << (i + 1));
        iFftButterflyNumPerStep = iFftSize / iFftPointNumPerBufferfly;
        for (j = 0; j < iFftButterflyNumPerStep; j++)
        {
            iTempStartIdx = j * iFftPointNumPerBufferfly;
            for (k = 0; k < (iFftPointNumPerBufferfly / 2); k++)
            {
                tCNTempWn = ptCNFftWn[k*iFftButterflyNumPerStep];
                tcn_mul(ptCNTempFftInputData[iTempStartIdx + k + (iFftPointNumPerBufferfly / 2)], tCNTempWn, &tCNTempMultData);
                tcn_plus(tCNTempMultData, ptCNTempFftInputData[iTempStartIdx + k], &ptCNTempFftOutputData[iTempStartIdx + k]);
                tcn_sub(ptCNTempFftInputData[iTempStartIdx + k], tCNTempMultData, &ptCNTempFftOutputData[iTempStartIdx + k + (iFftPointNumPerBufferfly / 2)]);
            } /* end of for for (int k = 0; k < (iFftPointNumPerBufferfly / 2); k++) */
        } /* end of for for (int j = 0; j < iFftButterflyNumPerStep; j++) */
        ptCNTempPointer = ptCNTempFftInputData;
        ptCNTempFftInputData = ptCNTempFftOutputData;
        ptCNTempFftOutputData = ptCNTempPointer;
    } /* end of for for (int i = 0; i < iFftSteps; i++) */
    for (i = 0; i < iFftSize; i++)
    {
        ptCNOutputData[i] = ptCNTempFftInputData[i];
    } /* end of for for (int i = 0; i < iFftSize; i++) */
    delete[]ptCNTempFftInputData;
    //free(ptCNTempFftInputData);
    delete[]ptCNTempFftOutputData;
    //free(ptCNTempFftOutputData);

    return 0;
}

/******************************************************************************
* 模块名称：ifft_butterfly(TComplexNumber ptCNInputData[], TComplexNumber ptCNOutputData[], TComplexNumber ptCNFftWn[], int piReshapeIdx[], int iFftSize, int iFftSteps)
* 模块功能：iFFT算法
* 参数说明：ptCNInputData[]  - TComplexNumber复数类型蝶形运算第一级输入数据序列buffer指针
*		    ptCNOutputData[] - TComplexNumber复数类型蝶形运算最后一级输出数据序列buffer指针
*		    ptCNFftWn[]      - TComplexNumber复数类型蝶形运算需要使用的输入旋转因子序列buffer指针
*		    piReshapeIdx[]   - int类型蝶形运算需要使用的输入数据排序索引序列buffer指针
*		    iFftSize         - int类型输入FFT总点数
*		    iFftSize         - int类型输入FFT蝶形运算级数
* 返回值  ：执行成功标识，返回0表示执行成功
* 建立时间：2020-08-12
* 创建人员：陈文棵
* 修改日期：
* 修改人员：
*******************************************************************************/
static int ifft_butterfly(TComplexNumber ptCNInputData[], TComplexNumber ptCNOutputData[], TComplexNumber ptCNFftWn[], int piReshapeIdx[], int iFftSize, int iFftSteps)
{
    int i;
    TComplexNumber *ptCNTempFftInputData;
    ptCNTempFftInputData = new TComplexNumber[iFftSize];
    //ptCNTempFftInputData = malloc(iFftSize * sizeof(TComplexNumber));

    for (i = 0; i < iFftSize; i++)
    {
        ptCNTempFftInputData[i].dReal = ptCNInputData[i].dReal;
        ptCNTempFftInputData[i].dImag = -1.0 * ptCNInputData[i].dImag;
    } /* end of for for (int i = 0; i < N_fft; i++) */
    fft_butterfly(ptCNTempFftInputData, ptCNOutputData, ptCNFftWn, piReshapeIdx, iFftSize, iFftSteps);      //调用FFT蝶形运算
    for (i = 0; i < iFftSize; i++)
    {
        ptCNOutputData[i].dReal = ptCNOutputData[i].dReal / iFftSize;
        ptCNOutputData[i].dImag = -1.0 * ptCNOutputData[i].dImag / iFftSize;
    } /* end of for for (int i = 0; i < iFftSize; i++) */
    delete[]ptCNTempFftInputData;
    //free(ptCNTempFftInputData);

    return 0;
}

template<class T>
int AMCUtilityAlgFFT::real_fft_power_dB(T pdInputData[], int iInputDataLength, T pdOutPutFftPower[], int *iOutputBufferLength, int iWantedFftSize)
{
    int iDvalue = 0;
    int iFftSteps = 0;
    int iFftSize = 0;
    int iError = 0;
    int i;

    /* 将输入的实数数据扩展为复数数据，数据不够的补0，用于FFT蝶形运算输入 */
    TComplexNumber *ptCNFftInputData;
    /* 计算FFT蝶形运算需要的输入数据重新排序索引 */
    int *piFftReshapeIdx;
    /* 计算FFT蝶形运算需要的输入旋转因子 */
    TComplexNumber *ptCNFftWn;

    /* 计算最接近的FFT蝶形运算级数及最接近数据点数的FFT点数 */
    iError = get_log2N(iWantedFftSize, &iFftSteps, &iDvalue);
    iFftSize = iWantedFftSize + iDvalue;


    //if (iFftSize > (*iOutputBufferLength))
    if ((iFftSize / 2) > (*iOutputBufferLength))
    {
        return -1;
    }/* end of for if (iFftSize > iOutputBufferLength) */

    ptCNFftInputData = new TComplexNumber[iFftSize];
    //ptCNFftInputData = malloc(iFftSize * sizeof(TComplexNumber));

    for (i = 0; i < iFftSize; i++)
    {
        if (i < iInputDataLength)
        {
            ptCNFftInputData[i].dReal = pdInputData[i];
            ptCNFftInputData[i].dImag = 0.0;
        }/* end of for if (i < iInputDataLength) */
        else
        {
            ptCNFftInputData[i].dReal = 0.0;
            ptCNFftInputData[i].dImag = 0.0;
        }/* end of for else */
    }/* end of for for (int i = 0; i < iFftSize; i++) */

    piFftReshapeIdx = new int[iFftSize];
    //piFftReshapeIdx = malloc(iFftSize * sizeof(int));
    iError = get_fftinput_reshape_idx(iFftSize, iFftSteps, piFftReshapeIdx);

    ptCNFftWn = new TComplexNumber[iFftSize / 2];
    //ptCNFftWn = malloc((iFftSize / 2) * sizeof(TComplexNumber));
    iError = get_fft_wn(iFftSize, ptCNFftWn);

    /* FFT蝶形运算 */
    iError = fft_butterfly(ptCNFftInputData, ptCNFftInputData, ptCNFftWn, piFftReshapeIdx, iFftSize, iFftSteps);

    /* FFT蝶形运算输出进行[0,2π)到[-π，π)的fftshift操作，并求样点功率 */
    for (i = 0; i < (iFftSize / 2); i++)
    {
        //pdOutPutFftPower[i] = ptCNFftInputData[i + iFftSize / 2].dReal * ptCNFftInputData[i + iFftSize / 2].dReal + ptCNFftInputData[i + iFftSize / 2].dImag * ptCNFftInputData[i + iFftSize / 2].dImag;
        pdOutPutFftPower[i] = ptCNFftInputData[i].dReal * ptCNFftInputData[i].dReal + ptCNFftInputData[i].dImag * ptCNFftInputData[i].dImag;
        if (pdOutPutFftPower[i] > 0.0)
        {
            pdOutPutFftPower[i] = 10 * log10(pdOutPutFftPower[i]);
        }
        else
        {
            pdOutPutFftPower[i] = -200;
        }
    }

    *iOutputBufferLength = iFftSize / 2; //将实际的FFT点数赋值给输出功率谱buffer长度，表示buffer中的有效点数

    delete[]ptCNFftInputData;
    delete[]piFftReshapeIdx;
    delete[]ptCNFftWn;

    return 0;
}


#endif // AMCUTILITYALGFFT_H
